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Abstract 

We present results for the phase diagram of three flavor QCD for //£<500 MeV. 
Our simulations are performed with imaginary chemical potential \x\ for which 
the fermion determinant is positive. Physical observables are then fitted by 
truncated Taylor series and continued to real chemical potential. We map out 
the location of the critical line T c ((j,b) with an accuracy up to terms of order 
(ub/T)®. We also give first results on a determination of the critical endpoint 
of the transition and its quark mass dependence. Our results for the endpoint 
differ significantly from those obtained by other methods, and we discuss possible 
reasons for this. 



1 Introduction 



The last two years have seen significant progress in simulating QCD at small baryon 
densities. Standard Monte Carlo methods fail in the presence of a non- vanishing chem- 
ical potential, for which the fermion determinant is complex and prohibits importance 
sampling with positive weights. While this problem remains unsolved for QCD, there 
are presently three avenues to investigate the (/i, T)-phase diagram, as long as the 
quark chemical potential in units of temperature, fi/T, is sufficiently small: (i) A two- 
dimensional generalization of the Glasgow method P predicts a phase diagram with a 
first order phase transition at large /i, terminating in a critical endpoint (ii) Taylor 
expanding the observables and the reweighting factor leads to coefficients expressed in 
local operators and thus permits the study of larger volumes [3]. (Hi) Simulations at 
imaginary chemical potential are not limited in volume since the fermion determinant 
is positive. They allow for fits of the full observables by truncated Taylor series, thus 
controlling the convergence of the latter, and subsequent analytic continuation to real 

The results for the location of the pseudo-critical line T c (/jlb) are consistent among 
all three approaches for baryon chemical potentials /iB up to /ab<500 MeV (hb = 3^i). 
This latter number gives the range of applicability of method (Hi), and hence the 
range over which convergence of the Taylor series can be checked explicitly. A review 
comparing these methods in detail can be found in [Sj. However, the location of the 
endpoint of the phase transition has only been computed by one method and for one 
set of quark masses [2] . Since the numerical study of critical phenomena is notoriously 
difficult even for fi = 0, it is important to cross-check this result by another method. 

In the present paper we investigate QCD with three degenerate quark flavors by 
lattice simulations with imaginary chemical potential. In this case we know reasonably 
well the chiral critical point m c (/x = 0), i.e. the critical bare quark mass m c for which 
the deconfinement transition changes from first order to crossover [HI 13 El- This point 
marks a second order phase transition in the universality class of the 3d Ising model 
[Tj- For the (/x,T) -phase diagram this means that for m < m c the line of first order 
deconfinement /chiral phase transitions extends all the way to the temperature axis at 
\i = 0, whereas for m > m c it terminates at a critical point (T*,fi*). The critical 
chemical potential, fi*(m) > 0, is expected to grow with the quark mass. Inverting 
this relation yields the change in the critical bare quark mass with chemical potential, 
m c (n). Our goal is to determine this function for chemical potentials /i£<500 MeV. 

We begin by computing the quark mass and chemical potential dependence of the 
critical line, T c (//B,m). Our simulations are accurate enough to allow for a determina- 
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tion of the (fi%) coefficient of its Taylor series. We then proceed to extract m c (fj,) by 
measuring the Binder cumulant of the chiral condensate. We find the /z-dependence 
of the latter to be very weak, proving a quantitative determination of m c (/i) to be 
extremely difficult. Just like the critical temperature, to c (/j) is an even function of 
the chemical potential with a Taylor expansion in fi 2 . We are only able to determine 
the first non-trivial coefficient with an error of 40%. However, this allows us to give a 
conservative upper bound on this coefficient at the 90% confidence level. This bound 
is in disagreement with a preliminary result obtained by Taylor expanded reweighting 

ESI- 

After summarizing the imaginary chemical potential approach in Sec. EJ Sec. 0] dis- 
cusses general qualitative features and expectations about the three flavor phase dia- 
gram and how it can be determined by simulations at imaginary chemical potential. 
Sec. |U introduces the Binder cumulant and its finite volume scaling as our computa- 
tional tool to determine the order of the phase transition as a function of quark mass 
and chemical potential. Our numerical results are presented in Sec. followed by a 
discussion and comparison with other work in Sec. |H1 Finally, we give our conclusions. 



2 The QCD phase transition from imaginary fi 

This section serves to fix the notation and summarize what is needed in the sequel. 
For a detailed discussion of the formalism we refer to [3]. The QCD grand canonical 
partition function Z(V, /J, T) = Tr ^e ~S/^<2/)/ T ^ ? w jth the same chemical potential 
/j for all flavors, can be considered for complex chemical potential /i = /ir + i/i/. Two 
general symmetry properties can be used to constrain the phase structure of the theory 
as a function of fij: (i) Z is an even function of fi, Z(p) = Z(—ft), where \x = fi/T; 
(ii) A non-periodic gauge transformation, which rotates the Polyakov loop by a center 
element but leaves Z unchanged, is equivalent to a shift in /i/ 0: 

Z(p R ,p, I ) = Z(p, R ,p, I + 27r/N). (1) 

For QCD (N = 3), these two properties lead to Z(3) transitions at critical values of 
the imaginary chemical potential, /if = %E (n + separating regions of parameter 
space where the Polyakov loop angle ((p) falls in different Z{2>) sectors. The resulting 
phase diagram in the (/i/,T) plane is periodic and symmetric about flj, as depicted 
in Fig. 1 for Nf = 3. The Z(3) transitions are first order for high temperatures, 
terminating at some critical temperature T c which coincides with the deconfinement 
line. This endpoint also belongs to the 3d Ising universality class [3]. The curves in 



2 







1/3 

iy(*r) 



2/3 



Figure 1: Schematic phase diagram in the (/xj,T) plane. Solid lines mark first or- 
der phase transitions, dotted ones crossover. The vertical line corresponds to a Z(3) 
transition and the curves to the deconfinement/ chiral transition at imaginary fi. Both 
terminate in critical points, each belonging to the 3d Ising universality class. 

Fig. 1 correspond to the deconfinement line continued to imaginary chemical potential, 
T c (fjji). We will map out these curves as well as the location of the critical endpoint 
as a function of quark mass. Their relation to real fi is discussed in detail in the next 
section. 

On the lattice, phase transitions can be studied by measuring fluctuations of gauge 
invariant operators, 



where we use the plaquette, the chiral condensate and the modulus of the Polyakov 
loop for 0(x). A transition region is signalled by a peak in the susceptibilities 



whose maximum implicitly defines the critical parameters, Xmax = xi^Pc)- In a 
finite volume, the susceptibility is always an analytic function of the parameters of 
the theory, even in the presence of a phase transition. The latter reveals itself by a 
divergence of Xmax in the infinite volume limit, whereas Xmax stays finite in the case 
of a crossover. This fact was used in [3] to show that, on any large but finite volume, 
the pseudo-critical coupling can be represented by a Taylor series in /x 2 . Here we also 
wish to study the quark mass dependence, and thus consider an additional expansion 
in m about the \x = chiral critical point m c (0). Hence we have for the pseudo-critical 






X = VN t ((SO) 2 ), 



(3) 



3 



coupling on a finite volume 



(3 c (afx, am) = ^2 c ki (a^) 2k (am — am c (0)) 1 . 

k,l=0 



(4) 



With the help of the two-loop lattice beta-function, the critical coupling can be con- 
verted to the critical temperature T c (m,/j,). In our previous work [I] we demonstrated 
by simulations that the //-series converges fast and the critical coupling viz. temper- 
ature are well described by the leading /x 2 -term. Analytic continuation between real 
and imaginary chemical potential is then trivial. 

Following [El HUE!, our strategy thus consists of measuring observables at (/2# = 
0,/2/ 7^ 0), fitting them by a Taylor series in Jx\ and then continuing the truncated 
Taylor series to real Jx. The Z(3)-transition closest to the origin, at pf\ = |, defines the 
convergence radius of the expansion and limits the prospects of analytic continuation. 
In physical units this corresponds to /x# ^,500 MeV. 

3 Qualitative features of the phase diagram 

In principle, as proposed in the order of the phase transition, and hence the location 
of the critical point, can be determined from a finite size scaling analysis of the critical 
coupling f3 c {V) itself, which attains its infinite volume limit as 



where a = 1 for a first order phase transition, a — 1/dv < 1 for a second order phase 
transition, and a = for a crossover. The critical endpoint on the curve is then defined 
by (3* = /3 c (/2*), for which the transition is of second order. However, such an analysis 
is not practical for analytic continuation. In a Taylor expansion in fl, the volume 
dependence resides in the coefficients, 



Clearly, a /i-dependent finite volume behavior as in Eq. (jSJ) cannot in general be well 
approximated by only a few terms of this series. 

Fortunately, the problem can be approached differently by considering variations of 
the quark mass, as outlined in |12j . In this case we have a three-dimensional parame- 
ter space {T,fi,m}. The critical temperature T c (/i, m) now describes a surface in this 
space, and the critical endpoint traces out a line T*{m) = T c (fi* (m) , m) , or equiva- 
lently T*(/x) = T c (m c (/x), /z), on this surface. Projections of this situation onto the 




(5) 



&(//, V) = &(0, V) + Cl {V)f + c 2 (lV + . . . 



(6) 
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(T,fi), (T, m) and (m, /i)-planes are shown schematically in Fig. |21 and El The bottom 
line in Fig. El (left) corresponds to the situation depicted in Fig.[U for some quark mass 
m < m c (0). With increasing quark mass, the critical endpoint of the deconfmement 
line in Fig. Amoves towards \x = 0, which it hits for m = m c (0). On the other hand, 
for decreasing quark mass it moves to larger ///, until it meets the Z(3) endpoint at 
some mass rh. For m < rh the deconfinement and Z(3) transition lines are connected. 

This feature appears in Fig.|3]as the intersection of m c (/x) with the vertical Z(3)-line, 
which can be shown as follows. In the chiral limit, the chiral condensate represents 
a true order parameter which is strictly zero in the deconfmed phase, and there must 
be a true phase transition for all values fi 2 < 0. Thus the line separating first order 
from crossover, m c (/x), cannot hit the negative /i 2 -axis unless it has an unexpected 
non-analyticity there, implying the existence of some rh > 0. 

Since the (pseudo-) critical line is analytic, so is the line of endpoints T*(/i), and 
by elimination of T the same holds for m c (/i). These are again smooth functions with 
analytic continuations to imaginary /i, which one may hope to describe well in terms of 
only a few coefficients. In our practical calculation we attempt to map out the phase 
diagram Fig. El by computing the coefficients of 



Preliminary results for the leading coefficient as determined from Taylor expanded 
reweighting, have been reported in ^3], and we will discuss this result in comparison 
to ours in Sec. |H1 

3.1 Theoretical expectations 

Before continuing to describe our calculational tools, let us make a few remarks about 
what one would expect theoretically. The analytic continuation approach has by now 
been tested for screening masses in the plasma phase [11] as well as for T c (/x) IT%] . 
In ^T] it was remarked that the screening masses are most "naturally" expanded in 
(n/ (vrT)), where "natural" means that the coefficients in such a series are of order one. 
The same observation is made regarding the critical temperature for the two flavor 
case. The result quoted in |2j can be rewritten as 



where the coefficient is of order one. In thermal perturbation theory this is easy to 
understand, as in the imaginary time formalism one expands in terms of Matsubara 



am c (/i 2 ) = ^c n (ajj) 



2n 



(7) 
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Figure 2: Left: critical lines in the (T ', /i 2 ) -plane for different quark masses m. The 
bold curve T*(/x) characterizes second-order transitions, separating the crossover and 
the first order regimes. Right: critical lines in the (T,m) -plane for different chemical 
potentials /i 2 . The bold curve represents T*{m). 




Figure 3: Schematic line of critical quark mass separating the first order and crossover 
region. The line is constrained by the \i = data point (diamond, Q/j and the fact that 
for m = the phase transition has to be first order for all imaginary fi 2 < 0, implying 
that intersection with the Z (3) -line happens at some quark mass m > 0. 

modes and the chemical potential always appears in this combination ^HE|- It is also 
transparent non-perturbatively in the case of an imaginary chemical potential /i/: the 
chemical potential gives an extra factor exp(i/x//T) for the boundary condition on the 
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fermionic fields, so that it is equivalent to shifting the Matsubara frequencies (2k+l)7iT 
by fij. Hence the relevant expansion parameter is the relative shift ((j,i/(ttT)). We will 
empirically confirm this for the three flavor case, where we also measure the next-to- 
leading coefficient. 

The same considerations apply to the quark mass expansion and provide a reason 
for the near independence of T c (n) upon the light quark masses [Sj. Fermionic modes 
contribute with non-zero Matsubara frequencies, and light quark masses are always 
negligible compared to those modes, (m/irT) <C 1. This is still the case for the strange 
quark mass, and so we expect the curve T c (/i) to be approximately the same for Nf = 3 
and N f = 2 + 1. 

For the critical quark mass one then similarly expects to have 



with ci of order one. The bare quark mass is not a physical quantity, but depends 
on the lattice action. For example, comparing calculations with p4-improved and 
unimproved staggered fermion actions, one finds m c (0)\ impr ps 0.25m c (0)\ unimpr [Zj. 
However, the mass renormalization should not be affected by // ^ 0, which is just an 
external thermodynamic parameter without ultraviolet renormalization. Multiplicative 
mass renormalization should therefore cancel out in the ratio Eq. Q, which, up to 
additive corrections 0(a 2 ), is directly comparable between different lattice actions. 

A remarkable finding for the critical temperature is that it is quite accurately de- 
scribed by the leading fi 2 term, at least up to |/2| = /if, where for imaginary /a the Z(3) 
transition occurs. The same was found for screening masses ^T] and recently also for 
the pressure |15j . We thus expect similar behavior for m c (fi), which will be confirmed 
by our simulations. If m c (//) is well described by the leading term, its intersection 
with the Z(3)-line at a quark mass value m > 0, as described in the previous section, 

(c \ 2 
ttt) ) — 0> 

or c\ < 9. 

4 Cumulant ratio and finite volume scaling 

In order to find the boundary between the first order and crossover regime along the 
critical line, we use the Binder cumulant ^Hj of the chiral condensate, 





5 4 



((<W) 2 ) 2 ' 



(10) 
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In the infinite volume limit this quantity assumes a universal value at a critical point. 
In particular, this observable was used in [7] to locate the chiral critical point at /x = 0, 
am c (0) = 0.0331(12) for staggered fermions on an iV t = 4 lattice, and to identify its 
universality class as that of the 3d Ising model, for which B 4 rs 1.604. On a finite 
volume, this value receives corrections. It also receives corrections away from the 
critical point, which are positive for crossover and negative for first order behavior. 
Cumulants calculated on increasing lattice sizes for different parameters will intersect 
at some pseudo-critical value of the parameters, with the -B 4 -value at the intersection 
point converging towards its universal value. 

In order to explicitly assess the quark mass and /i-dependence, we fit our data by a 
Taylor expansion about the \x = critical point, 

B 4 (m, fjt) = b m ( am ~ am c (0)) n (an) 21 , (11) 

n,l 

with boo(V — > oo) = 1.604. This observable can also be directly related to the critical 
line in the phase diagram Fig. El At the expansion point m c (0) we have B A = boo, and 
this value is maintained along the line m c (/i), which is a line of constant B A . This line 
is implicitly defined by the equation 

B 4 (m c (fx), fjt) = b 00 , (12) 

and in particular one obtains the coefficient c\ of Eq. © through the chain rule 

dam c _ dB A f dB A \ _1 _ _V 
<i(a/i) 2 d(afi) 2 \dam J bio 

If the volumes are large enough, the approach to the thermodynamic limit is governed 
by universality. In this case the volume dependence hidden in the coefficients of the 
series can be made explicit. Approaching the critical endpoint, the correlation length 
diverges as £ ~ r~ v ', where r is the distance to the critical point in the plane of 
temperature and magnetic field-like variables, and v = 0.63 for the Ising universality 
class. In practice, we first find (3 C for a given pair (m, /i), and then compute _B 4 for 
those values of the couplings. Since (3 = (3 C always, we thus have r = \m — m c (/i)|. 
£> 4 is a function of the dimensionless ratio L/£, or equivalently (L/^) 1 ^. Hence one 
expects the scaling behavior 

B, ((L/0 1/u ) = B A (L x ' v {am - am c (/x)) . (14) 
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5 Numerical results 



In our simulations we consider QCD with the Wilson gauge action and three degenerate 
flavors of staggered fermions, with bare quark masses in the range 0.025 < am < 0.04. 
We monitor finite volume scaling behavior using three lattice sizes, 8 3 x 4, 10 3 x 4 and 
12 3 x 4. The Monte Carlo employs the R- algorithm JK] with a step size St = 0.02, 
which is sufficiently small for the systematic errors 0(5r 2 ) to be negligible compared to 
our statistical errors. For each simulated parameter set we accumulate 10k — 70k unit- 
length trajectories, measuring the gauge action and the Polyakov loop and estimating 
the first four powers of the chiral condensate after each trajectory. The pseudo-critical 
values P c (afj,i) are obtained from a range of typically 4 simulated /3-values by means 
of the Ferrenberg-Swendsen reweighting method Hence, every data point in the 
following figures for the critical coupling and the cumulant ratio typically consists of 
over 100k trajectories. 

5.1 The critical line T c (/i) 

The calculation of the critical line proceeds as in the two flavor case [3J. The critical 
coupling j3 c was determined by finding a peak in the plaquette susceptibility, and we 
have checked that the chiral condensate and the Polyakov loop give consistent results. 
Our first task then is to determine the coefficients in the Taylor expansion Eq. (J3J). 

In Table^we give an exhaustive list of all possible three, four and six parameter fits to 
our data. For the expansion point in the quark mass, we have chosen am c (0) = 0.0323, 
which will be the result obtained in Sec. 15.31 Apart from resolving the leading linear 
quark mass and quadratic chemical potential dependence, our statistics is now also 
large enough to permit some statements concerning the next-to-leading terms. On our 
L = 10 lattice we studied the largest a/xj, and consequently get the most constrained 
fits for the /x 4 -term. Note that on this volume a quartic term is required to fit the 
data, while the other possibilities give significantly worse fits. The situation on the 
other volumes is consistent with this. The best four parameter fit is in all cases the 
one with a /i 4 -term. The other options give coefficients that are either consistent with 
zero within 1.5 standard deviations, or inconsistent between the different volumes. On 
the other hand, six parameter fits do not significantly reduce the x 2 ; and thus are not 
fully constrained yet. 

Comparing the coefficients of the fits including the quartic term between the volumes, 
we observe that the present statistics is unable to resolve systematic finite volume ef- 
fects, all volumes being compatible within one standard deviation. It is then expedient 
to further constrain the fit parameters by fitting all volumes together. We use the best 
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L 


coo = &(0, 0) 


Cio, (At 2 ) 


C20, (/U 4 ) 


On, (m) 


c 2, ijn 2 ) 


c u , (/i 2 m) 


rVdof 


8 


5.1450(2) 


0.790(8) 




1.748(26) 




_ 


0.94 


8 


5.1451(2) 


0.739(26) 


0.93(46) 


1.735(27) 





_ 


0.66 


8 


5.1449(3) 


0.791(8) 





1.745(27) 


2.3(6.3) 


_ 


1.01 


8 


5.1450(2) 


0.789(8) 





1.751(31) 





-0.085(0.40) 


1.02 


8 


5.1449(4) 


0.740(27) 


0.99(48) 


1.712(44) 


6.8(9.9) 


0.38(0.63) 


0.75 


10 


5.1449(2) 


0.800(5) 


_ 


1.831(23) 


_ 


_ 


7.66 


10 


5.1457(2) 


0.667(15) 


1.98(22) 


1.818(23) 





_ 


1.21 


10 


5.1457(2) 


0.798(5) 





1.808(23) 


-28.2(4.6) 


_ 


5.41 


10 


5.1449(2) 


0.798(5) 





1.883(40) 





-0.38(25) 


8.04 


10 


5.1458(2) 


0.679(17) 


1.78(25) 


1.820(42) 


-8.6(5.6) 


-0.07(26) 


1.21 


12 


5.1455(4) 


0.764(13) 


— 


1.770(44) 


— 


— 


1.04 


12 


5.1457(5) 


0.721(31) 


0.94(63) 


1.788(46) 


— 


— 


0.44 


12 


5.1454(5) 


0.763(13) 


- 


1.784(57) 


7.5(18.9) 


- 


1.48 


12 


5.1443(8) 


0.791(21) 




2.28(31) 




-2.5(1.5) 


0.12 


8-12 


5.1442(1) 


0.7943(33) 




1.796(16) 






4.10 


8-12 


5.1453(1) 


0.705(10) 


1.46(15) 


1.780(16) 






1.43 


8-12 


5.1453(2) 


0.7917(34) 




1.792(16) 


-14.2(3.3) 




3.68 


8-12 


5.1448(1) 


0.7940(34) 




1.807(23) 




-0.1(16) 


4.20 


8-12 


5.1457(2) 


0.705(10) 


1.43(16) 


1.767(26) 


-10.1(3.9) 


0.10(19) 


1.17 



Table 1: Fits of the Taylor expansion /3 c (m, /i) ; Eq. Q) ; to our data. 



four parameter fit highlighted in the table as our final result for the critical coupling, 
which is shown in Fig. H]as a function of // 2 . 

We conclude that we have a signal for a /i 4 contribution to the critical coupling. This 
is a result of having more accurate data and does not invalidate our earlier observation 
that the line is well described by the leading term. E.g. at fij the contribution of this 
term to the critical coupling is only ~ 0.1%. Converting to continuum units by means 
of the two- loop beta-function as in [3], we thus obtain for the critical line in three flavor 
QCD 

T c (/i,m) n „, /17 > fm-m c (0Y 



1 + 1.937(17) 



T\T, 



-0.602(9) 



+ 0.23(9) 



(15) 



T c (0,m c (0)) 

T c on the right side of this and the following two equations is meant to be the same 
as in the denominator on the left, which plays a role for the coefficient of next-to- 
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5.14 




P c -c 01 (am-am c (0)) 
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5.08 







0.01 0.02 0.03 0.04 0.05 0.06 0.07 



Figure 4: Combined data for L = 8 — 12 and various quark masses. Data for different 
am are shifted to am c (0) according to the best fit of Table 1, which is also shown. A 
weak /i 4 -dependence is visible. 

leading terms. Note that the mixing of the errors on the parameters in the critical 
coupling through standard error propagation drowns out the /i 4 -signal in continuum 
units. Since we do not yet have a signal for a mixed (m/i 2 )-term, the quark mass and 
chemical potential dependence are separately consistent with 



A mixed dependence only appears in higher orders, having no effect at our present 
accuracy. This explains why critical lines obtained previously for various different 
quark masses agree so well [S]. 

We may then directly compare our result with existing ones for Nf = 2 jl] and 
Nf = 4 [Hj in Fig. EJ As one would expect, our result falls between these two. Note, 
however, that these earlier results were not sensitive to a /z 4 -term, which makes itself 
felt at the right end of the interval and also lowers the /^-coefficient, cf. Tabled Also 
shown in the figure is the result for Nf = 2 + 1 as obtained by reweig hting 2J. In 
accordance with our expectations from Sec. 13.11 due to the quark mass independence 
of the critical line this result is practically identical to the one for Nf = 3. 



T c (/x,m c (0)) 

Teh*, m ) 
T c (0,m) 



T c (/i,m) 
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0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 

JJ. E /T liJT 

~B c r b c 

Figure 5: Left: One sigma error bands on T c (/ie) for different Nf (Nf = 4 from 
Only the Nf = 3 calculation is accurate enough to include a quartic term. Right: 
Comparison of Nf = 3 with Nf = 2 + 1 from J^j. 



5.2 First order vs. crossover and error estimates 

Before presenting our results for the cumulant ratio, we make some remarks concerning 
the considerable technical difficulty of these measurements. Inspection of the Monte 
Carlo history of an observable over a sufficiently long Monte Carlo time reveals that 
the tunneling frequency between the different vacua is very low: observing only one 
crossing per a few thousand trajectories is typical. This is expected in a first order 
regime, where tunneling is suppressed by a potential barrier, but the same observation 
is made in the crossover regime. 

The reason for this behavior is the fact that, on the lattice sizes used here, the 
probability distributions for measurements at the critical coupling (3 c (m, ii) have not 
yet reached their asymptotic scaling regime. This is illustrated in Fig. |H1 where we 
show the distributions of plaquette values on two volumes for a point each in the first 
order and crossover regimes. In accordance with expectation, the first order region 
displays a two peak structure and tunneling gets more suppressed on a larger volume. 
In the crossover region we observe accordingly a merging of the two peak structure with 
increasing volume. However, this merging to the asymptotic Gaussian distribution is 
not yet complete, and a remnant of the two-peak structure can be clearly identified. 
The displayed parameter values are deep in the crossover region, and the situation gets 
only worse closer to the critical point. This is a well known difficulty in the investigation 
of phase transitions. 
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Figure 6: Distribution of plaquette values at (3 C for L — 8, 10 and /i = 0. Left: First 
order transition, am = 0.025. Right: Crossover, am = 0.04. 

On the other hand, the value of B^ and its statistical error are driven by the number 
of tunnelings rather than the total number of measurements. Essentially the observable 
distinguishes between crossover and a first order transition by picking up the difference 
in the frequency of tunnelings. This leads to a much slower reduction of error bars than 
in the case of the critical couplings, where only the change of the observable between 
the two phases is needed, for which the number of measurements is relevant. Hence, too 
short Monte Carlo runs with less than a few tens of tunnelings tend to underestimate 
the statistical error on B4. More dangerous is the finite volume remnant of tunneling 
suppression in the crossover regime which can, for too short runs and combined with 
too small an error estimate, lead to an underestimate of B4 and hence to misidentifying 
a crossover as a weakly first order signal. 

In light of this, we can only be fully confident of our B4 error estimates for L = 8 
lattices, where tunneling is faster and we have the longest Monte Carlo runs. On this 
volume we obtain a significant result for the /i-dependence, whereas for L = 10, 12 the 
signal is hidden in the noise. These volumes will be mainly used for consistency and 
scaling checks. 

5.3 The cumulant ratio as function of m and fi 

Following Sec. El we proceed to discuss our measurements of the cumulant ratio B4 
along the critical line in order to determine its endpoint and its quark mass dependence. 
Our current accuracy constrains only the leading terms 0(am, (a/i) 2 ) in the Taylor 
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am am 



Figure 7: Binder cumulant for (3 = (3 C as function of quark mass for afii = (left) and 
a/ij = 0.2 (right). 

expansion of B4. To begin, we perform an analysis analogous to the one in For 
a fixed value of a/xj, we measure _B 4 along the critical line p c (am, a/Xj), cf. Fig. |2 
(right). The critical quark mass separating first order from crossover behavior is then 
extracted from the intersection of B4 measured on different volumes. This is shown 
in Fig. 0for a\ij = and a/xj = 0.2. Our results for a/x/ = are in full agreement 
with those reported in serving as a check of the analysis. The volume dependence 
appears to be moderate, and for the intersection point between the larger volumes we 
get am c (0) ~ 0.033, compared to am c (0) = 0.0331(12) However, practically the 
same result is obtained for a\ii = 0.2, pointing to a very weak //-dependence of B±. 
Indeed, plotting our data for fixed am as a function of (a///) 2 , no structure beyond 
noise is apparent to the eye. 

In order to obtain better accuracy we modify our analysis. Let us rewrite the leading 
terms of the Taylor expansion Eq. (jllj) as 

B^am, a/x) = 1.604 + B (am — am c (0) — A(a/x) 2 ) , (17) 

where we have traded the parameters {b 00 , 6 i, &10} f° r { m c{fy, A, B}. With the con- 
stant fixed to its infinite volume value, finite volume corrections to 600 wm now show 
up in m c (0), which can be compared with the previous result. In this form we can 
collapse all our data obtained for various pairings (am, a///) into one plot and fit them 
by a single three parameter fit. Finally, d(am c ) / d(afi) 2 , as in Eq. (|13p. is now immedi- 
ately given by the fit parameter A. Fig. |H] shows all L = 8 data combined in this way 
together with the best fit. The fit results for all volumes are displayed in Tabled 

However, even after combining all data on one volume, the /^-coefficient A is still 
only weakly constrained. The data on the larger lattices are consistent with a negligible 
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T 1 =1 




0.025 0.03 0.035 0.04 



(am+ACa^j) ) 

Figure 8: Binder cumulant in the (am, a/ii) -plane for L = 8, the line represents the fit 
given in Tabled 
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0.044(19) 
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0.0325(3) 
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0.54 
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0.0320(5) 


-0.010(10) 


25.4(2.1) 


0.53 


12 


0.0326(2) 




35.8(1.7) 


0.19 


12 


0.0325(4) 


-0.008(18) 


35.0(2.5) 


0.22 


16 


0.0331(3) 




57.0(6.3) 


0.19 


L 


am c (0) 


A 


B/L 1 '" 


X 2 /dof 


8-12 


0.0323(3) 


0.0008(8) 


0.67(3) 


0.84 



Table 2: Fits of the Taylor expansion B4, Eq. ([771 ), to the data. L = 16 data for \i = 
are taken from |7J/. 



yU-dependence, as is apparent by the acceptable fits obtained without such a term. Only 
on the 8 3 lattice, for which we have the best statistics, do the fits prefer a positive value 
of this quantity. 

Let us now try to combine the different volumes by exploiting the fact that m c (0) 
appears close to its infinite volume limit, and hence -B 4 should be close to the scaling 
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Figure 9: Finite volume scaling of the fit parameter B from Eq. The line repre- 

sents a fit to ~ L x l v , with v = 0.62(3), x7 dof = °- 2 - 




-0.2 0.2 0.4 0.6 
(am-am c (0)+A (au/)L (1/v) 

Figure 10: Combined fit for all volumes by fixing the scaling to the Ising form, Eq. 

region on the volumes considered. In order to explicitly test for this, we plot the fit 
parameter B against the volumes for which it was obtained, and fit the data to the 
expected asymptotic scaling behavior B 4 (L) ~ L 1 ^, cf. Eq. (JT3J). For this purpose, we 
also use the L = 16 data from [7j. This is shown in Fig.El and the resulting v = 0.62(3) 
is indeed consistent with the Ising value v ~ 0.63. Having thus established the explicit 
volume dependence of B^ as in Eq. (|14j) . we may combine all available volumes into 
one maximally constrained fit, in order to get higher accuracy. This is done in Fig. [TU1 
The resulting parameter values are given in the last line of Table El Our result for the 
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zero density critical quark masss then is am c (0) = 0.0323(3), with a statistical error of 
1%, and in perfect agreement with the calculations at /i = [7J |K| • 



5.4 The line of critical endpoints 

The situation is less clear for the //-dependence, as Table |2] shows. The combined fit 
over all volumes does not constrain the parameter A enough to yield a non-vanishing 
result. Clearly, this is due to the L = 10, 12 lattices, whose results are consistent with 
zero, but whose negative central values neutralize the significant answer obtained on 
L = 8. Since these lattices only add noise to the determination of A, we thus quote 
the L = 8 number as our tentative final result, 



with higher terms being smaller than our present error of 40%. A check of the fit result 
is obtained by measuring for different quark masses along the vertical Z(3)-line 
at afij, cf. Fig. 01 in order to determine m. We have done so on L = 10 and find 
0.029 <am< 0.032, while Eq. (JUJ) in lattice units predicts am = 0.030. 

Note that, in terms of our natural expansion units, the coefficient of interest is 
not unnaturally small. Determining it to better accuracy is, however, a formidable 
numerical task that requires computational resources on the largest scales available. A 
more conservative result is obtained by adding two standard deviations to the central 
value, resulting in a bound c\ < 1.6 at 90% confidence level. 

Taken at face value, Eq. (fT%|) tells us the critical bare quark mass for a given chemical 
potential as sketched in Fig. El while its inverse yields the location of the critical 
endpoint for a given bare quark mass. The renormalization of the bare quark mass 
cancels in the ratio, so that it should be independent of the lattice action chosen, 
up to additive cut-off effects. Moreover, since in the mass range of interest the zero 
temperature pion mass M% oc m, we have 



where M^(0) « 290(190) MeV for unimproved (p4-improved (preliminary)) staggered 
fermions, respectively . These numbers highlight the strong need to eliminate 
cut-off effects on M=(0). 

Another result involving only physical quantities is obtained by eliminating the bare 
quark mass in computing the line of critical endpoints, 




(18) 



{M C M? = (M^(0)) 2 



m c {n) 



(19) 



m c (/i = 0) ' 




(20) 
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While in this function we lose the information on the quark mass dependence, the 
curve relates only infrared quantities and describes a physical property of the QCD 
parameter space, cf. Fig. |21 

6 Discussion 

Let us now try to compare our results to those obtained by other groups. The first 
Taylor coefficient in m c (/i 2 ) for the three flavor theory was also calculated by means of 
Taylor expanded reweighting ^H] • In the form of Eq. (fTH|) their result for the coefficient 
is 67(19), compared to ours of 0.82(36). We observe that, continued to imaginary 
this result violates the bound c\ < 9 derived in Sec. 13.11 The only way to avoid this 
conclusion would be large (9(/i 4 )-effects, for which we see no evidence. While at present 
we have no explanation for this rather drastic disagreement, we speculate that it is a 
statistics problem: the preliminary result of ^3] is based on six thousand trajectories, 
and measurements for different \i are always correlated in reweighting approaches. 
Considering the problems we mentioned in Sec. 15.21 the similarly sobering findings of 
Ref . |B] , and the scatter of our uncorrelated data in Fig. |H1 this might account for the 
discrepancy. 

Eventually, we are of course interested in the 2+1 flavor theory with non-degenerate 
masses. In this case the line of constant B4 derived from the leading order expression 
Eq. (dU) reads 

2(m M - m c (0)) + (m s - m c (0)) - V = . (21) 

However, a linear extrapolation in the quark mass to am s is most likely not valid. 
Blindly substituting the bare quark masses of j2] and our value for A, one would obtain 
a critical chemical potential fis ~ 3 GeV. While this number is certainly meaningless, 
it seems nevertheless that our calculation would put the critical endpoint of the de- 
confinement line at considerably larger values of \ib than those reported in EI] • To 
avoid extrapolations over large ranges, 2+1 flavor QCD with physical masses requires 
additional calculations in the light and heavier quark mass regimes and could be quite 
different numerically. 

Finally we would like to add one more comment concerning the difficulties of distin- 
guishing a first order phase transition from crossover, Sec. 15.21 Our discussion focused 
on the observable B4, and one may ask about its relevance for other methods of de- 
termining the endpoint, like finite size scaling of susceptibilities or Lee- Yang zeroes. 
While other observables might well have smaller statistical errors than B± when mea- 
sured on the same number of configurations, their relative behavior between first order 
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and crossover regimes is nevertheless driven by the number of tunnelings, and there- 
fore suffers from the same slowness of the Monte Carlo history as our B 4 measurement, 
requiring similar statistics in order to arrive at reliable results. 

7 Conclusions 

We have investigated the finite density phase diagram of three flavor QCD for yU#<500 
MeV by means of lattice simulations at imaginary chemical potential. Compared to 
previous studies with Nf = 2, 4, we gathered much increased statistics allowing us to 
determine the location of the critical line T c (m, [Ib) through terms linear in the quark 
mass and quartic in the chemical potential. The curvature of the critical line becomes 
more negative with increasing Nf. Any mixing terms between quark mass and chemical 
potential are smaller than our present accuracy, rendering T c (m, [iB)/T c (m, 0) quark 
mass independent to a good approximation. 

We have also studied the nature of the phase transition along the critical line, and the 
location of its endpoint as a function of quark mass, by studying the Binder cumulant 
as a function of quark mass and chemical potential. We were able to compute the first 
coefficient of the critical quark mass m c (/z 2 ) to 40% accuracy. A constraint at the 90% 
confidence level puts our result at considerable odds with a preliminary result given by a 
Taylor expanded reweighting technique [15J, our critical endpoint being at larger /j,b for 
comparable quark masses. Our central results are given in Eqs. (|T5 |) .([Tfi )) and (fT%j) . (f2Tty . 
While we have clearly demonstrated the feasibility of such a calculation, our results 
exhibit the formidable difficulty of this task, whose unambiguous completion requires 
computational resources beyond those presently available to us. An extrapolation to 
the physical 2 + 1 flavor case requires additional simulations to account for the heavier 
strange quark, and is envisaged for the future. 
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